### This script plots the analyses to compare different ways to incorporate the PDG
### into the Zonation analyses

### Libraries
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.5.2
library(dplyr)
## Warning: package 'dplyr' was built under R version 3.5.2
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)
## Warning: package 'tidyr' was built under R version 3.5.2
### Read data

## Read output data generated by the script spatial_analyses_zonationVSproxiesdivgen_all_WT_final.R
## this shows the area per species and proxies of genetic diversity within it 
## considering the Zonation solution for 20% of the country considering only 
## species distribution models (no habitat, etc).
sols_summary_spp<-read.csv("../data/comparations_output/sols_summary_spp.txt")

# check data
head(sols_summary_spp)
##                Solution Richness  shannon   simpson   Area Prop_to_AreaSP
## 1                ENM_sp       64 3.216956 0.9341559 473887      1.0000000
## 2          Zonation_MDP       58 2.908383 0.9179912 127904      0.2699040
## 3       Zonation_MDP_ZV       59 2.929987 0.9210733 127082      0.2681694
## 4      Zonation_MDP_PDG       61 3.565616 0.9677778 142830      0.3014010
## 5   Zonation_MDP_vs_PDG       64 3.718850 0.9720398 124221      0.2621321
## 6 Zonation_MDP_PDG_ADMU       63 3.182794 0.9418937 129665      0.2736201
##   mean.prop median.prop Dist.to.Optimun                            sp
## 1        NA          NA               0 Capsicum annuum glabriusculum
## 2 0.3786268   0.2580767              NA Capsicum annuum glabriusculum
## 3 0.4083096   0.3092418              NA Capsicum annuum glabriusculum
## 4 0.5480085   0.4495986              NA Capsicum annuum glabriusculum
## 5 0.6525528   0.6421699              NA Capsicum annuum glabriusculum
## 6 0.4407066   0.3329768              NA Capsicum annuum glabriusculum
### Analyses and plots

## Check mean proportion of species area and mean proportion of proxies within it, per solution
group_by(sols_summary_spp, Solution) %>% 
  summarise(., mean.prop.sp.area =  mean(Prop_to_AreaSP) , mean.prop.proxies = mean(mean.prop))
## # A tibble: 6 x 3
##   Solution              mean.prop.sp.area mean.prop.proxies
##   <fct>                             <dbl>             <dbl>
## 1 ENM_sp                            1                NA    
## 2 Zonation_MDP                      0.476             0.543
## 3 Zonation_MDP_PDG                  0.367             0.576
## 4 Zonation_MDP_PDG_ADMU             0.481             0.663
## 5 Zonation_MDP_vs_PDG               0.408             0.757
## 6 Zonation_MDP_ZV                   0.429             0.555
### Plots

## Violin plots

# dist to optimun

filter(sols_summary_spp, Solution!="ENM_sp") %>%
  ggplot(., aes(x=Solution, y=Dist.to.Optimun, color=Solution)) + geom_violin() + geom_jitter(size = 0.3) +
  stat_summary(fun.y=mean, geom="point", color="red") +
  theme(axis.text.x = element_blank()) +
  scale_y_continuous(name="Distance to Optimun")
## Warning: Removed 595 rows containing non-finite values (stat_ydensity).
## Warning: Removed 595 rows containing non-finite values (stat_summary).
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf
## Warning in min(diff(sort(x))): no non-missing arguments to min; returning
## Inf
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf
## Warning in stats::runif(length(x), -amount, amount): NAs produced
## Warning: Removed 595 rows containing missing values (geom_point).

# Simpson

filter(sols_summary_spp, Solution!="ENM_sp") %>%
  ggplot(., aes(x=Solution, y=simpson, color=Solution)) + geom_violin() + geom_jitter(size = 0.3) +
  stat_summary(fun.y=mean, geom="point", color="red") +
  theme(axis.text.x = element_blank()) +
  scale_y_continuous(name="Simpson Index")

# Area sp conserved

filter(sols_summary_spp, Solution!="ENM_sp") %>%
  ggplot(., aes(x=Solution, y=Prop_to_AreaSP, color=Solution)) + geom_violin() + geom_jitter(size = 0.3) +
  stat_summary(fun.y=mean, geom="point", color="red") +
  theme(axis.text.x = element_blank()) +
  scale_y_continuous(name="Proportion of distribution of the sp conserved")

# Mean prop area proxies by sp ponserved

plot_b <-filter(sols_summary_spp, Solution!="ENM_sp") %>%
  ggplot(., aes(x=Solution, y=mean.prop, color=Solution)) + geom_violin() + geom_jitter(size = 0.3) +
  stat_summary(fun.y=mean, geom="point", color="red") +
  theme(axis.text.x = element_blank()) +
  scale_y_continuous(name="Mean proportion of the area of each proxy conserved by sp")

plot_b

# Median prop area proxies by sp ponserved

filter(sols_summary_spp, Solution!="ENM_sp") %>%
  ggplot(., aes(x=Solution, y=median.prop, color=Solution)) + geom_violin() + geom_jitter(size = 0.3) +
  stat_summary(fun.y=mean, geom="point", color="red") +
  theme(axis.text.x = element_blank()) +
  scale_y_continuous(name="Median proportion of the area of each proxy conserved by sp")

# Median and mean together

filter(sols_summary_spp, Solution!="ENM_sp") %>%
  gather(., SummaryStat, value, -Solution, -c(Richness:Prop_to_AreaSP), -Dist.to.Optimun, -sp) %>%
  ggplot(., aes(x=interaction(SummaryStat, Solution), y=value, color=SummaryStat)) + geom_violin() + geom_jitter(size = 0.1) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Proportion of the distribution of the sp conserved vs MEAN proportion of each proxy conserved by sp

filter(sols_summary_spp, Solution!="ENM_sp") %>%
  ggplot(., aes(x=Prop_to_AreaSP, y=mean.prop, color=Solution)) + geom_point() +
  scale_x_continuous(name="Proportion of distribution of the sp conserved") +
  scale_y_continuous(name="Mean proportion of the area of each proxy conserved by sp")

# Proportion of the distribution of the sp conserved vs MEADIAN proportion of each proxy conserved by sp

filter(sols_summary_spp, Solution!="ENM_sp") %>%
  ggplot(., aes(x=Prop_to_AreaSP, y=median.prop, color=Solution)) + geom_point() +
  scale_x_continuous(name="Proportion of distribution of the sp conserved") +
  scale_y_continuous(name="Median proportion of the area of each proxy conserved by sp")

## with regression lines

plot_a <-filter(sols_summary_spp, Solution!="ENM_sp") %>%
  # plot points
  ggplot(., aes(x=Prop_to_AreaSP, y=mean.prop, color=Solution)) + geom_point(size=0.7) +
  scale_x_continuous(name="Proportion of species distribution (%)") +
  scale_y_continuous(name="Mean proportion of area of proxy represented by species distribution (%)", 
                     breaks = seq(0, 1, by = 0.25), expand = c(0, 0)) +
  # plot fitted curve
  geom_smooth(method=loess, aes(fill=Solution)) 
plot_a

### Plot for paper figure:


plot_a <- plot_a +
  # nicer background
  theme_bw()+ theme(panel.border = element_blank(), panel.grid.major = element_blank(),
                    panel.grid.minor = element_blank(), 
                    axis.line = element_line(colour = "black")) +
  
  # nicer legend (fill and color are needed becasue we are plotting both points smooth and points)
  scale_fill_discrete(name = "Scenarios",
                      labels = c("SDM (n=116)", 
                                 "SDM+LZ (n=143)", "SDM+PGD (n=218)",
                                 "SDM*PGD (n=5004)", "SDM and PGD as ADMU (n=117)")) +
  scale_color_discrete(name = "Scenarios",
                       labels = c("SDM (n=116)", 
                                  "SDM+LZ (n=143)", "SDM+PGD (n=218)",
                                  "SDM*PGD (n=5004)", "SDM and PGD as ADMU (n=117)")) +
  theme(legend.position = c(.99, 0.01), legend.justification = c("right", "bottom")) +
  
  # title
  labs(title="a)") +

  # larger text
  theme(
  plot.title = element_text(size=14, face="bold"),
  axis.title = element_text(size=14),
  axis.text = element_text(size=13),
  legend.text = element_text(size=13),
  legend.title = element_text(size=13, face="bold"))

plot_a

plot_b <- plot_b +
  # nicer background
  theme_bw()+ theme(panel.border = element_blank(), panel.grid.major = element_blank(),
                    panel.grid.minor = element_blank(), 
                    axis.line = element_line(colour = "black")) +
  # no legend 
  theme(legend.position = "none") +
  
  # nicar x names
  scale_x_discrete(name="Scenarios",
                   labels=c("SDM", 
                            "SDM+LZ", "SDM+PGD",
                            "SDM*PGD", "as ADMU")) +
  # title
  labs(title="b)") +
  
  # larger text
  theme(
    plot.title = element_text(size=14, face="bold"),
    axis.title = element_text(size=14),
    axis.text = element_text(size=13))

plot_b  

## Plot Together


# Multiple plot function taken from http://www.cookbook-r.com/Graphs/Multiple_graphs_on_one_page_(ggplot2)/
#
# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)
# - cols:   Number of columns in layout
# - layout: A matrix specifying the layout. If present, 'cols' is ignored.
#
# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
# then plot 1 will go in the upper left, 2 will go in the upper right, and
# 3 will go all the way across the bottom.
#
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)
  
  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)
  
  numPlots = length(plots)
  
  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  
  if (numPlots==1) {
    print(plots[[1]])
    
  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    
    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

multiplot(plot_a, plot_b, cols=2)